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Abstract 

RNA viruses are known to replicate with extremely high mutation rates. 
These rates are actually close to the so-called error threshold. This threshold 
is in fact a critical point beyond which genetic information is lost through 
a second-order phase transition, which has been dubbed the "error catastro- 
phe." Here we explore this phenomenon using a field theory approximation 
to the spatially extended Swetina-Schuster quasispecies model [J. Swetina 
and P. Schuster, Biophys. Chem. 16, 329 (1982)], a single-sharp-peak land- 
scape. In analogy with standard absorbing-state phase transitions, we develop 
a reaction-diffusion model whose discrete rules mimic the Swetina-Schuster 
model. The field theory representation of the reaction-diffusion system is 
constructed. The proposed field theory belongs to the same universality class 
than a conserved reaction-diffusion model previously proposed [F. van Wi- 
jland et al, Physica A 251, 179 (1998)]. From the field theory, we obtain 
the full set of exponents that characterize the critical behavior at the error 
threshold. Our results present the error catastrophe from a new point of view 
and suggest that spatial degrees of freedom can modify several mean field 
predictions previously considered, leading to the definition of characteristic 
exponents that could be experimentally measurable. 

PACS numbers: 87.10. +e, 02.50.-r, 64.60.Ak 



1 



Typeset using REVTgX 



I. INTRODUCTION 



RNA viruses offer a unique opportunity for exploring long term evolution under con- 
trolled conditions due to their high mutation rates |]||. Their evolutionary success is to a 
large extent due to their small-sized genomes, but specially to their enormous plasticity and 
adaptability to changing environments ||. These viruses display the highest possible muta- 
tion rates and, as a consequence, their populations, the so-called molecular quasispecies [|J, 
are extremely heterogeneous. The quasispecies structure has numerous implications for the 
biology of viruses. The most relevant of them is that mutant swarms are reservoirs of vari- 
ants with potentially useful phenotypes in the face of environmental change. As extremely 
simplified entities at the border of a life-like state, they are specially apt to mathematical 
modeling |)]J§. 

The replication of these molecules implies two basic reactions MM: (i) error- free copy, 
when a (molecular) species Jj replicates by using available monomers (A), i.e.: 

( A ) + I .±Si2I i , (1) 

and (ii) mutation: 

(A) + 1^ It + Ij (i^j). (2) 

The parameters Ai and Qi are the replication rate and the quality factor, respectively. 
Qi G [0, 1] is a measure of the correctness of the replication process, and it is maximum 
(Qi = 1) if no mutations occur. ^ are the mutation rates which can lead to transitions 
between species j — > i. 

The standard approach to the quasispecies dynamics (in the limit of very large pop- 
ulations) is based on the continuous Eigen model ||. Here a set of molecules which can 
replicate and mutate is considered. The basic equations are: 

^ = (AiQi - Di) Xi + ^ijX 3 + $ i5 (3) 

where Xi, with i = 1,2, ...,n, accounts for the population size of each species, stands 
for spontaneous degradation of molecules (assumed to be linear), and $j is an outflow term 
which takes into account the removal of molecules from the system. If we introduce the 
constraint of constant population size (CP), J2i x i — const., the previous equations read: 

^ = (AiQi - A - E) Xi + ]T VijXj, (4) 

where the mean value of the so-called excess productivity, — Ai — Qi, is given by E = 
Tl,i(Ai — D{)xi/ J2j x j- 

One of the most important results of Eigen's theory was the finding that a phase tran- 
sition takes place when mutation rates are tuned. Specifically, let us assume for simplicity 
that each species is composed by a string of elements {Si, S u } of size v M. We consider 
in total a population of iV strings. A particular kind of sequences, composed by the ele- 
ments {Si , ...,S^}, represent the correct genomic sequences of the species, and is called 
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the master sequence. Each time a string is chosen to be replicated, it does so with some 
sequence-dependent probability r, = P({Si}). If replication occurs, each unit can mutate 
with probability /i, and reproduce exactly with probability 1 — /i. Mutation introduces 
disorder into the system, and it can be shown that a well-defined mapping exists between 
replication dynamics and the two-dimensional Ising model in such a way that the tempera- 
ture is given by T « — | log(/i/(l — fi))\ P-pH- For small /i the replication system reaches 



a steady state in which the probability of observing the master sequence is finite. On the 
other hand, for values of \x larger than a given threshold, the probability of observing the 
master sequence is vanishing small. 

Eigen's theory predicts that genetic information becomes lost for mutation rates higher 
than the critical rate /i c , due to a breakdown of heredity and the lack of selection — the so- 
called "error catastrophe." It has been shown that this phenomenon indeed occurs in RNA 
viruses, which replicate close to the error threshold 0. Actually, experimental data reveal 
that real viruses have mutation rates /i ~ 1/z/, that is, inversely proportional to the size of 
the genomic content, consistently with the prediction, and this has led to the claim that 
increased mutation rates might be able to bring virus populations into extinction. Such a 
strategy has been recently shown to hold in vitro and is likely to be feasible in vivo f!2 |. 



The presence of a critical mutation rate allows to interpret the error catastrophe in the 
framework of standard absorbing-state phase transitions (APT) fI5,H|. APT are a class of 



non-equilibrium transitions in which, by the variation of a control parameter, the system 
crosses from an active phase with everlasting activity, to an absorbing phase, in which the 
system remains trapped forever, with no possibility to escape. In the framework of species 
replication with mutation, the active phase is identified with the low mutation regime, 
while the absorbing phase corresponds to the high mutation regime. Most APT are phase 
transitions of second order. If we characterize the system by an appropriate order parameter 
t/>, which in this case corresponds to the density of master sequences in the system, by tuning 
the parameter /i we observe the typical behavior 

ip = 0, for /i > fi c , 

ip ~ Oc - /J,) 13 , for n < fj, c , 

close to the critical point fi c . The previous expression serves to define the critical exponent j3. 
The analogy with error catastrophe is in this sense clear: for mutation rates larger than the 
error threshold, the virus is inviable and it quickly dies. For small mutation rates the virus is 
able to survive and reaches viable populations whose size is an increasing function of \l c — /i. 
Further extending the analogy with APT, we can consider the spatial and time dependence of 
the order parameter if), and define the correlation function g(r, t) = (ip(r' , t')ip(r' + r,t' + t)), 
where the bracket denote averages over different realizations of the system. According to the 



dynamic scaling ansatz ||15|| , we expect to observe close to the error threshold the behavior 

(5) 



which defines the correlation length £, related to the distance to the /i c by 

^(/i c -/i)-^. (6) 
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Eqs. (|5|) and define the new critical exponents rj, v± and z, which determine the scaling 
of the correlation function with respect to changes in the mutation rate \i. 

Based in the previous analogy, in this paper we propose to study the phenomenon of 
the error catastrophe from the point of view of an APT by analyzing a reaction diffusion 
model with captures the essence of the replication-plus-mutation mechanism of quasispecies 
dynamics. The model allows the construction of an associated field theory, representative 
of the same universality class, following a standard technique outlined in the work of Doi, 
Cardy, and others [|IT| -|T9|j . The field theory developed here is shown to correspond to a 
conserved reaction diffusion model previously proposed by van Wijland et al |21J (see also 
Ref. ||21[|), in which the critical exponents were obtained performing a one- loop renormaliza- 
tion group analysis. 

An important aspect, seldom considered in previous studies, is the effect of spatial degrees 
of freedom in quasispecies models. An exception is Adami's work on artificial life systems, in 
which a set of replicating bit strings of code spread on a two-dimensional lattice [2^] . Under 
appropriate conditions, it was shown that the population spontaneously evolves to the error 
threshold, although no characterization of the model behavior at this critical point was 
performed. Besides, only a few experimental studies have recently reported the presence of 
several patterns of virus distribution that cannot be explained in terms of spatially-implicit 
quasispecies dynamics H23|j24} . 

The interest in this problem is twofold: On the one hand, virus populations show hetero- 
geneity in space, thus introducing further complexity in quasispecies dynamics and creating 
new opportunities to viral evolution. On the other hand, it would be important to know if 
spatially extended, mean-field models, are appropriate descriptions of the real quasispecies 



dynamics in space. Although RNA fitness landscapes are known to be rugged [£5|,[2(J, here 
we consider the simplest, single-sharp-peak landscape |27]. This model has been used as a 
null model of quasispecies populations and a field theory of the reaction-diffusion rules can 
be developed. 

The paper is structured as follows. In Sec. II we propose, in analogy with the Swetina- 
Schuster model |27].|28| a streamlined react ion- diffusion model which captures the minimal 
elements in the reproduction/ mutation mechanism of quasispecies dynamics. Sec. Ill re- 
ports a mean-field analysis, which allows to pinpoint the key parameters of the model. In 
Sec. IV we construct the field theory corresponding to the model. In analogy with the 



analysis performed by van Wijland et al. and Kree et al. [p0,21|, we obtain the relevant 



critical exponents. Finally, we interpret our results and put forward several experimental 
applications of them in Sec. V. 



II. REACTION-DIFFUSION MODEL 

The key ingredient of our model consists in the assumption that one of the sequences 
B = I m has a high replication rate, while all the others A = Ij^ m have the same, lower 
replication rate. The first sequence is called the master sequence and this approximation 



defines the so-called Swetina-Schuster model P7| . Second, by assuming that the sequences 



are long enough (consistently with real RNA viruses, where v > 10 4 ), backward mutations 
from A to B can be neglected p9[ . 
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In order to propose the reaction steps defining the model, we consider a simplified version 



of the Swetina-Schuster model (see Ref. [28 1 and references therein). In this model, a popu- 
lation of N strings of size v evolves by a mechanism of replication with errors. Each string is 
defined by a sequence S%, . . . S u , with Si = {0, 1}. At each time step, we select a string and 
replicate it, after removing another string chosen at random. The replication takes place 
with probability 1 for the master sequence, defined by Si = 1, Vi, and with probability p < 1 
for the rest. The replication procedure replaces each element Si of the string for S[ = Si 
with probability 1 — fi, and for S^ = (Si + 1) mod 2 with probability /x. This version of 
the model shows in plain view the elementary steps of the error catastrophe: replication of 
the master sequence at a certain rate, mutation of the master sequence, and lack of back- 
ward mutation, for sufficiently large sequence length. Another important ingredient in the 
model to be remarked is the implicit constraint of constant number of sequences, realized in 
the random deletion step, and which is usually implemented in the quasispecies analytical 
models Q . The model described in Ref. seems to display all the features of the error 
catastrophe. 

Using this simple framework, we can translate the dynamics of the Swain-Schuster model 
at a microscopic level in terms of reactions among particles of type B and A, corresponding 



to the master sequence and the mutants, respectively. The simplified model in Ref. 
implicitly introduces interactions among sequences, by means of the random deletion step. 
Thus, our reaction diffusion model considers all the possible binary reactions between par- 
ticles of type B and A, that are compatible with the outcome of the rules used in |28) . The 
set of reactions that we consider is: 

B + B B + A (7a) 

B + A B + B (7b) 

A + B A + A (7c) 

The steps represent the replication/mutation of the first species, coupled with the random 
deletion of the second species. Thus, the reaction ( |7a| ) implements the replication with mu- 
tation of a master sequence B, which happens with an effective mutation rate fi, coupled 
to the random deletion of a sequence of type B\ the reaction ( fH5| ) represent the exact repli- 
cation of a master sequence, at rate 1 — /i, with the deletion of a A sequence; finally, the 
reaction (173) stands for the exact replication of a sequence A, at rate A, together with the 
deletion of a master sequence. All the remaining binary reactions with replication/mutation 
plus random deletion do not alter the total number of particles, and are thus not consid- 
ered. The proposed set of reactions mimic the conserved nature of the model imposed by 



the random deletion of sequences in Ref. J28[, in a more nature way than in the original 



quasispecies model, Eq. (|3[), in which one had to impose an external flow term in order 
to ensure conservation. 

The set of equations (0) constitutes the core of our model. Spatial effects are taken into 
account by allowing the different particles to diffuse with respective diffusivities Da and Db 



Given the interpretation of the different particles, it is natural to consider Da = Db, 
that is, both master sequence and mutants diffuse with the same speed. However, for the 
sake of completeness, we will develop the formalism with Da ^ Db, and make them equal 
only as a last step. 
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III. MEAN-FIELD ANALYSIS 



In order to gain some preliminary intuition on the behavior of the model, we analyze 
it by applying a standard mean-field analysis. Let us denote by pB and pa the density of 
species B and A, respectively. Since the reactions (0) conserve the number of particles, the 
total density ps + Pa is constant in time: 

Pb + Pa = P- (8) 
The classic (mean-field) equations for the densities ps and pa are readily found to be: 

?£jL = (1 — A - p)pbPa - ppl, (9a) 

rtn * 

-(1 - A - p)pbPa + PP B - (9b) 



dt 

Combining this with the conservation condition (|8|) we obtain a single equation for the 
density ps of master sequences: 

^ = (l-X-p)p PB -(l-X)p 2 B . (10) 

This equation has two stable stationary states, depending on the value of p: 

p B = 0, for /i > 1 - A, (11) 
Pb= l \ X _~/ p for/i<l-A. (12) 

At the mean-field level we observe the presence of a standard absorbing-state phase transition 
at a critical point p c = 1 — A. In the subcritical regime, p > p c , the order parameter (in 
this case the density of master sequences) vanishes; in the supercritical region, p < p c , the 
order parameter has a power-law dependence on p: 

p B ~ (p c - pf, (13) 

which defines the critical exponent (3 in the mean- field approximation, /3mf = 1- As we will 
see in the next section, the presence of fluctuations will change the value of (3 at the relevant, 
experimental, dimensions. A very interesting property of this conserved RD model is that 
the critical point is independent of the total particle density p, and is given as a function 
only of the reproduction rate A. This situation should be compared with the conservative 
RD systems proposed so far, in which the total particle density plays the role of the tuning 
parameter, and must be tuned to a critical density p c in order for the system to display 
critical behavior pp] , pT . 



The mean-field solution also provides the expression for the probability V that there is 



at least one master sequence in the steady-state regime |32fl : 

V = Q(p c -p), (14) 

where O is the Heaviside function. 
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IV. FIELD THEORY 



Our model consists in a set of particles of type B and A diffusing in a hypercubic lattice 
of mesh size h and size L, that interact probabilistically according to the rules (|7|) whenever 
they meet at the same lattice site. The dynamics of the model is defined through a master 
equation for the probability P({n},{m},£) of having a particle configuration {n}, {m} 
of particles B and A, respectively, at time t. The configuration {n} = {n!,n 2 , . . . ,n L d} 
represents the occupation number of each node in the lattice. The master equation for P is: 



^P({n},{m},t) = 

%r E tK- + m 



rii - 1, rij + 1, 



{m},t)-n l P({n},{m},t)] 



+ ^ E [K + !) p (W> . . . m< - 1, m 3 - + 1, . . . , t) - m i P({n}, {m}, *)] 



(15a) 

(15b) 
(15c) 



+ A* E + 1 ) n i- p (- • • ^ + 1, . . . mj - 1, . . . , i) - - l)P({n}, {m}, £)] 

i 

+ (1 - n) [(rii - 1)K + 1)P(. . . , m - 1, . . . , mi + 1, . . . , t) - n imi P({n}, {m}, t)} (15d) 

i 

+ ^ E + 1 )( m « _ ^-^C- • • , + 1, . . . , mj - 1, . . . , t) - nimiP({n}, {m}, t)] , (15e) 



where Db and are the diffusion coefficients for B and A particles, % is summed over all 
the lattice sites, and j over the nearest neighbors of the site i. The first two terms in the 
rhs of Eq. (|15|) implement diffusion through a random hopping of particles between nearest 
neighbor sites. The initial condition P({n}, {m},0) is given by a Poisson distribution, with 
an average density per site equal for both types of particles. 

The next step consists in recasting the master equation into a "second quantized" form, 
following the procedure described by Doi ITT6HT9T . We introduce two sets of annihilation and 
creation operators at each lattice site, Oj and b\ for B particles, and a, and a\ for ^4 particles, 
which fulfill the standard commutation rules 



Sij. 



(16) 



With this commutations rules, the operators have a bosonic character, that is natural given 
the multiple occupancy of sites allowed in the model. With the help of the vacuum state |0), 
defined by hi |0) = 6, |0) = 0, we construct an orthonormal basis of states \n,m), defined by 



n, m 



(17) 



and work in the Fock space spanned by this basis. In terms of this Fock space, the state of 
the system at time t is represented by the vector state \P(t)), defined as 



l^)>= E P({n},{m},t)\n,m). 

{n},{m} 



(18) 
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In terms of this vector state, the master equation Eq. fll5|) can be rewritten as a Schrodinger 
equation in imaginary time 



0_ 

Of 



\P(t)) = -H\P(t)) 



with a Hamiltonian, or time-evolution operator, H defined by: 

h = E ~ a}) (a* - + ^(SJ - &})& - 6, 



(19) 

(20) 
(21) 



Eq. 



can be formally solved in terms of the operator H yielding 

\P(t)) = exp(-Ht) |P(0)) . 



(22) 



From this solution, it is possible to derive all the statistical properties of the RD system, 
applying a projection technique [16 19| . For practical purposes, it is convenient to map this 
second-quantized form into a field theory, using a coherent state representation. Performing 
a time-slicing of the evolution operator in Eq. fl22"|), via the Trotter formula, we can express 
the vector state \P(t)) as a path integral, weighted with the exponential of an action S, over 
a set of classical fields a*, a, b*, and b, which are related with the two types of particles. After 
taking the continuum limit (h — > 0), the vector state can be written as the path integral 
over space and time dependent fields 

\p(t)) = Jva Va* Vb Vb* exp(-%, a*, b, b*}) |P(0)) , (23) 

where the action S has the form [|33| 

S[a, a*, 6, b*} = J d d x J dt [a*[d t - D A V 2 }a + b*[d t - D B V 2 ]b 

+ fi{b* - a*)b*b 2 + (1 - fi){a* - b*)b*ab + \{b* - a*)a*ab) . 

Within this formalism, we can compute the average value of any observable P({n}, {m}) 
performing the path integral 



(F(t)} = C Jva Va* Vb Vb* F(a, b) exp(-S[a, a*, 6, &*]), 

where C is an appropriate normalization constant. 

The final step in the derivation of the field theory consists in performing the shift 



(24) 



1 + 5, 



and the change of variables 



t = b, 
ip = b — a, 



1 + 6, 



a + b — p, 
a. 



(25) 



(26) 
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The final action describing the RD system is 



J d d x J dt {ij[d t ij - D B V 2 tfj -r^ + g^ 2 - g 2 ^} + 0[<9 t - D A V 2 (j> + 7 V 2 ^] (27) 
+ ^[-g^ - viip(j) + v 2 ip 2 } + tp4>[-g 4 ip - v 3 ip(j) + v A ip 2 }} , 
where we have defined the coupling constants 



r 


= 9i = 


(1 -X-n)p, 


9i 


= Vi = 


1-A, 


92 


= v 3 = 


1 - A - n, 


93 


= pv\ = 


= (i-m)p> 


v 2 


= 1, 
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= D A - 


-D B . 



The coupling constant r, 7, and Vi are coarse-grained versions of the microscopic reaction 
rates. Since we are only interested in the behavior of the system close to a critical point, 
however, the actual value of these parameters is irrelevant [ B^ET] ]. 

A naive power counting shows that the critical dimension of this field theory is d c = 4, 
and that the coupling constants i>j are irrelevant, and can be in principle discarded. With 
this final form, it is easy to recognize that the action (|27D represents the same field theory 
analyzed by van Wijland et al. for the conserved reaction diffusion model 

B + A — >2B, 
B — > A. 



In their study, the authors worked out the renormalization group analysis for this system, 
for both cases: Da < Db and Da = Db, providing the critical exponents up to a one-loop 
expansion. The case Da = Db, which also corresponds to the universality class of a model 
of population dynamics with pollution described by Kree et al. ||21||, is the relevant one in 



the problem under consideration. Quoting the results of Refs. [pO| , |2T| , we have the critical 
exponents: 

13 = 1 " k (28a) 

v± = —r (28b) 

V = ~ (28c) 
z = 2, (28d) 

where e = 4 — d gives the dimensionality of the system. The results for u± and z are exact, 
derived field-theoretically by analysing the symmetries of the action (|2T|), and are thus valid 
for all dimensions. The values of f3 and rj, on the other hand, and expansions around e — 0, 
and thus they are expected to hold only for small values of e. 
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In view of the results (|28|), the relevant exponents at the physical dimension d = 3 are 

"± = 2/3 z = 2, 

(3 ~ 0.969 r? ~ -0.125. 1 ; 

The values of u± and z are exact, while those for j3 and r\ represent an approximation given 
by the replacement of the small parameter e for 1. In dimension d = 2 or less, the exponents 
(p8|) have to be taken with a grain of salt, due to the terms t>j in Eq. (p7|), which might 
become relevant at low dimensions [201. 



V. DISCUSSION 

The dynamical theory of molecular evolution developed by Eigen and Schuster reveals 
the presence of an intrinsic, sharp limit to molecular information carriers. This threshold 
is a generic feature of replicator systems involving reproduction and mutation. The Eigen- 
Schuster theory predicts that, under the effect of evolutionary pressures selecting for high 
variability, such replicators will evolve towards the error threshold. This is the case of RNA 
viruses and experimental evidence clearly supports this theoretical prediction. 

Previous theoretical models have analysed the stochastic dynamics of quasispecies under 
different approaches. But all of them considered spatially- implicit models (mean-field-like), 
paying no attention to local effects derived from incomplete mixing. Here we have explored 
this problem using the simplest quasispecies model, described by a single-sharp-peak repli- 
cation landscape. The aim of our study was to see how the statistical behavior of a spatially 
extended molecular replication system would differ from the mean-field predictions. 

We have considered a simplified reaction-diffusion model where two types of "particles" 
(the master sequence and the mutant sequences) diffuse, replicate, and mutate on a given 
spatial domain. Applying the standard approach to absorbing-state phase transitions, a 
field theory has been developed and it has been shown to be the same reported by Wijland 
et al. |20| and Kree et al. [^T] . 



The main message from our study is that relevant differences between mean-field models 
and real dynamics are expected to be observed even in the simplest scenario considered 
here. Larger deviations should be expected in more realistic models incorporating a better 
description of molecular replicators and their dynamics. In particular, the sharpness of the 
transition (as defined by the f3 exponent) is not very different at different dimensions. This 
suggests that no measurable differences should be expected to be observed in experimental 
systems. The correlation exponent u±, however, does change appreciably, from u±(d = 1) = 
2 to Ux(d — 3) =2/3 thus leading a faster decorrelation at realistic dimensions. Such 
an increase in u± will enhance the coexistence of different strains (quasispecies) in a given 



spatial domain and thus the probabilities of success for the virus. 

Several caveats of this approach are worth mentioning. It is known that real RNA 
viruses have actually multipeaked landscapes P5] , pB| . We will expand our analysis to such 



situation in a future work (although the associated field theory, constructed in terms of 
several particles of type B^, each one representing a different viable master sequence, is 
expected to be much harder to develop and analyze). However, in many situations the 
quasispecies are observed to be confined in a fitness peak, so that our previous analysis 
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essentially holds |25|,|26[]. Also, the real RNA virus dynamics takes place through a virus- 
cell interaction not considered here, while several previous theoretical models used in order 
to understand well-defined experimental results (where a cell population was present) have 
been shown to be successful in providing a full understanding of the evolutionary dynamics 
of RNA populations |35|,|36 . 
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